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Abstract. - This contribution is devoted to a new model of HIV multiplication. We take into 
account the antigenic diversity through what we define "antigenicity" whether of the virus or 
of the adapted lymphocytes by a specific variable. We model two processes in the interaction 
of the immune system and the viral strains. On the one hand, the presence of a given viral 
quasi-species generates antigenically adapted lymphocytes. On the other hand, the lymphocytes 
kill viruses for which they have been designed. We consider also the mutation and multiplication 
of the virus. A new infection term is derived. 

So as to compare our system of differential equations with some well-knowns models, we 
study all of them mathematically and compare their predictions to ours in the reduced case of 
only one antigenicity. In this particular case, our model does not yield any major qualitative 
differences. We prove that in this case, our model is biologically consistent (positive fields) and 
has a unique continuous solution for long time evolution. In conclusion, this model improves the 
ability to simulate more advanced phases of the disease. 

Resume. - Nous proposons et etudions un nouveau modele pour la multiplication du virus du 
sida. Nous prenons en compte la diversite antigenique via la notion d' antigenicity tant pour 
les virus que pour les lymphocytes qui leur sont adaptes. Cette antigenicite est caracterisee 
par une variable specifique. Nous modelisons deux phenomenes d'interaction entre le systeme 
immunitaire et les souches virales. D'une part la presence d'une quasi-espece virale suscite des 
lymphocytes adaptes a leur antigenicite. D'autre part les lymphocytes tuent les virus pour 
lesquels ils ont ete produits. Nous prenons aussi en compte la mutation et la multiplication du 
virus. Nous proposons un nouveau terme traduisant l'infection. 

Afin de comparer notre systeme d'equations differentielles avec ceux existant dans la littera- 
ture, nous en faisons une analyse mathematique. Nous comparons leurs predictions avec la notre 
dans le cas simple d'une seule antigenicite. Dans ce cas particulier, notre modele ne presente pas 
de difference importante avec ceux de la litterature, est biologiquement consistant (les grandeurs 
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restent positives) et a une unique solution continue aux temps longs. En conclusion, ce modele 
ameliore la capacite de simulation des phases plus avancees de la maladie. 
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MSC Classification Codes: 34-99, 65-05, 65Z05, 92B99, 92C50 

1 Introduction 

Virus multiplication is at the basis of viral infection. Although the viral replication cycle 
involved makes heavy use of the infected cell's resources ( Warrilow et al, 2008; Fujii et al., 
2007), the enzyme (s) in charge of viral genome replication is (are) frequently encoded in the 
latter (Brass et al., 2008). Compared to cellular polymerases, viral polymerases are usually 
more error-prone (Kashkina et al, 2006). It follows that the viral mutation rate, defined 
as the average number of base changes at a given position of the genome per replication 
cycle, may be large compared to that observed in our own cells, which is about one in a 
billion (Petravic et al., 2008). For certain viruses, in particular retroviruses like HIV-1, 
it could be up to 1,000,000 times greater. By the end of 2007, the Los Alamos HIV data 
base listed over 230,000 different viral sequences (see www.hiv.lanl.gov). Given the small 
genome size of these viruses, chances are then that each member of a viral progeny carries 
mutations. A single point mutation may be enough to simultaneously affect two genes 
encoded in different, but overlapping reading frames, whilst silent mutations, which do 
not change amino acid coding, may nevertheless have important biological consequences 
(Mueller et al, 2006) 

A virus is mainly characterized by its ability to infect target cells (infectivity) and by 
its antigenic signature (antigenicity), defined as both the capacity to induce an immune 
response and also its strength and type. Immunogenicity is the ability of antigens to elicit 
a response from cells of the immune system. Mutations during virus replication may 
therefore release infective or non-infective viruses, of the same or of different antigenicity. 
For HIV-1, the ratio of infectious to non-infectious particles is estimated to range from 
1:1 to 1:60,000, depending on the type of cell infected and the viral strain (Thomas et 
al., 2007). Whether the virus is infective or not, over 800 mutations affecting HIV-1 
antigenicity were identified in its envelop gene (env) alone (Kothe et al, 2006). 

By encoding its own replication enzymes, the virus has control over its replication fi- 
delity and thereby challenges heavily the immune system, due to the huge burden imposed 
by the number of infective virions produced and their antigenic diversity. This burden 
is even worse when the virus targets part of the immune system (CD4 displaying cells) , 
as is the case for HIV-1. In addition, the immune cell proliferation induced by the viral 
attack will provide HIV-1 virions with new targets, engaging the cell-virus dynamics in 
an exponentially soaring extension regime. 

Kinetic modeling is therefore of high interest for understanding the course of infection. 
It is a prerequisite for designing and optimizing treatment strategies based on antiviral 
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drugs. A large number of deterministic and stochastic inter- and intra-cellular models of 
HIV dynamics have already been proposed (Finzi and Silliciano, 1998; Nowak and May, 
2000; Pastore and Zubelli, 2008; Perelson, 2002; Snedecor, 2003; Perelson and Nelson, 
1999; Tuckwell and Le Corfec, 1998), but none enable the prediction of the course of the 
disease in all its phases. For instance, even if, following antiviral treatment, the plasma 
load of the virus becomes undetectable, unscheduled bursts occur, probably fed by vi- 
ral sanctuaries disseminated in various tissues and organs (lymph or neuronal tissues, 
gastro-intestinal or uro-genital tracts etc. See Recher et al. (2007)). One may assume 
that, despite tissue-specific kinetic diversities, the course of infection in all sanctuaries 
(including plasma) obeys a common, complex host-predator relation, differing only by 
sanctuary-specific parameters. The course of the global infection would then be the result 
of all local processes. This result would however not be a simple addition or superposi- 
tion, as it is likely that each local process would provid e viruses having locally-specified 
antigenicities and infectivities which may challenge the immune cells in the same or other 
sanctuaries. 

What is the infection phenomenology ? As far as our study is concerned, the pro- 
cess involves four major participants or "fields": uninfected T lymphocytes (denoted T), 
infected ones (U), infectious viruses (V) and non-infectious viruses (U). For each par- 
ticipant, a characteristic antigenicity is recognized by a lymphocyte, or displayed by a 
virus. In the following, we call antigenicity the variable associated to this biomolecular 
characteristic and denote it by the index j. 

The evolution with time of Tj depends on 3 phenomena: regression of the Tj population 
due to the viral attack (whatever the antigenic pedigree of the infecting virus), which 
transforms T/s into U/s; the stimulation of the immune system by both infectious and 
noninfectious viruses of antigenicity j; the natural fate of Tj independently of the viral 
presence, ie spontaneous generation and death of Tj species. The Tj population may 
include cross-reacting species, which are active also against a (limited) number of targets 
with different antigenicities. 

Uj lymphocytes are derived only from the Tj population following viral attack, ie 
the immune system does not directly contribute inactive T cells of antigenicity j. It is 
assumed also, that the switch from the Tj to the Uj state occurs only upon viral infection. 
As for Tj, C/j's are subject to natural death. 

Viruses Vj and Wj are generated through infection of any susceptible T cell, what- 
ever its antigenicity. The parent virus of Vj and Wj may be a Vj (no viral antigenicity 
modification) or a V\. of different antigenicity (viral antigenicity mutation) . 

Viruses Vj and Wj will be the target of lymphocytes Tj exclusively. Both viral species 
have a natural death rate. 

Viruses Vj and Wj differ by mutations in genes involved in infectivity, but not affecting 
antigenicity j. 

This modelling assumes that viral genome parts responsible for infectivity may differ 
from those responsible for antigenicity. Since the viral strategy is to escape the immune re- 
sponse whilst minimizing loss of infectivity, excess mutations in the antigenicity-specifying 
part of the viral genome would be more favorable. The immune system senses mainly the 
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viral surface, hence mutations in the viral envelop genes would be most beneficial, but 
they should not, or marginally only, affect viral genome parts responsible for infectivity. 
HIV-1 handles in part this dilemma by introducing mutational hot-spots in its genome (Ji 
and Loeb, 1994 and the website [http:/ /www.hiv.lanl .gov), mainly in the envelop genes, 
where the mutation rate is much higher due to local sequence and structural particularities 
of the genome (Bebenek et al, 1993; Kothe et al., 2006). 

What is the therapeutic motivation ? 

Approved drugs for AIDS treatment are of three kinds mainly. Two of them inhibit 
reverse transcription (using nucleotide and non-nucleotide analogs) and the third inhibits a 
viral enzyme in charge of cleaving reverse transcriptase from a precursor protein. Because 
of its high mutation rate, HIV rapidly develops resistance to any one of these drugs taken 
individually. Resistance can be considerably delayed by using various combinations of 
these drugs (multitherapies). So far however, no combination has been found that could 
clear the virus. Therefore therapies are life-long and unfortunately have considerable 
side-effects. 

Obviously, the high mutation rate of reverse transcriptase is central to the successful 
viral strategy. It allows the virus to escape the circulating immune cells (antigenic muta- 
tions) and to develop drug resistance, although over 90% of its progeny lacks infectivity 
and will therefore be rapidly cleared. The natural viral mutation rate is at the limit of the 
"error threshold" (Biebricher and Eigen, 2005), as a slight increase would produce 100% 
non-infectious viruses. Conversely, reducing the mutation rate would reduce the antigenic 
diversity and allow the immune system to eliminate the stabilized viral strains, and drug 
resistance would vanish. 

New AIDS therapy at stake? 

A promising therapeutic approach would then be to take control of the viral mutation 
rate. This was shown to be feasible (Derrien, 1998 with a CNRS-filed patent based on 
this work USPTO 6,727,059; Murakami et al, 2005; Harris et al, 2005; Drosopoulos et 
al, 1998), by supplying the reverse transcriptase with nucleotide analogs. Some of them 
relax while others reinforce the replication fidelity, without blocking reverse transcription. 
For a review see Anderson et al. (2004). 

Both therapeutic strategies would give rise to specific viral dynamics. These need 
to be understood and assessed in detail. The medical decision to choose one particular 
strategy and setting up the adapted drug regimen for a patient, given his viral load and 
lymphocyte count (dose and extent of treatment, time expected to reach viral clearance 
etc) needs careful analysis. Simulation of viral dynamics with a drug regimen could help 
in reaching this decision. To this end, the present work is to build a realistic mathematical 
model of viral dynamics. 

In the following we review some well-known models by giving an analysis of the fixed 
points and their stability in Section 2. In Section 3, we derive our new model precisely 
and Section 4 is devoted to a full study of its mathematical properties in a reduced case 
(only one antigenicity). We conclude in Section 5. 
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2 Some popular models 

Throughout this section, we review some well-known models. Some of them take specific 
biological reality into account. So as to come to a common description with our model 
presented later, we reduced them in a preliminary step when needed. The fields T, U, V, W 
have the same meaning as above. When the model has a term identical to ours, we denote 
the parameter of the term as ours. When it is different, we denote it the same way as the 
authors and add a subscript depending on the authors. Also, we use the very same values 
of common parameters to have comparable results. We take most values from Snedecor 
(2003) since this paper is very well documented and check these values with other articles 
(Nelson et al, 2001; Herz et al, 1996; Perelson et al, 1996; ...). 

All the fields of the models are non-dimensionnalized by using the value of the non- 
infected lymphocytes at health (no virus and long time) as a characteristic value both for 
the lymphocytes (infected or not) and the virus. 

2.1 A Nowak-May model 

In Nowak and May (2000) the authors propose a simple model for the dynamics of the 
non-infected lymphocytes (T), the infected ones (U) and the viral population (V): 

' ^ = (5{l-T)- lNM VT, 

< ^ = lNMVT-aU, (1) 
[ ^ = aU-i NM V. 

This system is linear except for one interaction term. A simple analysis of this system 
gives two fixed points that are either health (T* = 1, U* = 0, V* = 0) or a state that we 
will call hereafter seropositivity since virus coexist with lymphocytes: 

T* = U* = ^(l - T*), V* = -^-(1 - T*). (2) 

Notice, that the link of our denomination with what is usually called "seropositivity" 
would need to be investigated. In our contribution, a seropositivity state is any fixed 
point different from health. We also define a seropositivity fixed point to be admissible if 
the fields (namely U*) are positive. Here it gives o^atm — o-Inm < and so T* < 1. We 
say that such a fixed point is immune deficient. So according to the parameters, there is 
either one fixed point (health) or two (health and seropositivity). 

Concerning stability, the study is easy too. The characteristic polynomial is — ((3 + 
A)(A 2 + (a + ^ata/)A + a^NM — o-Inm)- The second order polynomial has a positive 
discriminant. On the one hand, if at^NM ~ o-Inm < (there exists a seropositivy fixed 
point as T* < 1) then health is unstable in one direction and stable in the two others. On 
the other hand if oi^nm — o-Inm > 0, health is the only fixed point and is locally stable. 
We perform numerical simulations for which the following parameters were taken: 

(3 = 0.01 day" 1 , ^ NM = 0.0125 day" 1 , a = 0.7 day" 1 , a = 250 day" 1 . 



6 



Frangois Dubois, Herve VJ. Le Meur and Claude Reiss 



In FigurefU one may see the dynamics for the case of health alone (£nm = 10 day -1 ), 
and an initial value close to health (T = 1, Uq — Vo — 5%). As is proved by theory, 
health (alone) is stable. After about one day, the number of viruses is decreasing by one 
order of magnitude. During this period, the number of sound lymphocytes is changed 
only by 1 part over 10,000. It begins to rise again after about 8 days. 




Figure 1: Numerical results for Nowak-May model with three equations. Case where the 
health is stable. Short time evolution on the left (3 days) and phase plane for lymphocytes 
and viruses on the right for a 600 days evolution. 



By changing £^vm for £nm = 1- day -1 , we reach the regime where there are two fixed 
points. With an initial condition identical to above (T = 1, Uq = Vq — 5%), we numeri- 
cally check that seropositivity is stable as can be seen on Figure El First the number of 
viruses rises from 0.05 to more than 63, then it decreases to 0.53. During the same time, 
the number of healthy lymphocytes decreases to 0.035 typically then raises to 0.36. In the 
end, the dynamics oscillates around a fixed point that is the seropositive state (T ~ 0.25 
and V ~ 1 with our parameters). 

In Nowak and May (2000), the authors also proposed some other models taking the 
antigenicity into account. In Chapter 12 of their book, they propose various models, 
suggesting that "antigenic variation generates the long-term dynamics that give rise to 
the overall pattern of disease progression in HIV infection" (p. 124). Their "general idea 
was that the rapid genetic variation of the virus generates over time viral populations 
(quasispecies) which are more and more adapted to grow well in the microenvironment 
of a given patient" (p. 124). 

According to their best model, "The immune system and the virus population are in 
a defined steady-state only if the antigenic diversity of the virus population is below a 
certain threshold value. If the antigenic diversity exceeds this threshold [then] the virus 
population can no longer be controlled by the immunne system." (Nowak and May, 2000 
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Figure 2: Numerical results for Nowak-May model with three equations. Case where a 
seropositive state is stable. Short time evolution on the left (100 days) and phase plane 
for lymphocytes and virus on the right for a 600 days evolution. 



p. 125). Such a behavior seems meaningful. With our notations their model writes: 

= Vj(r NM - p NM Tj - q N MZ), Vi = 1, ...N 
= cnmVj - b NM T 3 -UV Tj, V* = 1, ...N 
= k NM V — bjqMZ — u NM V Z. 



d^ 

dt 
d7} 

dt 
dZ_ 

~dt 



where V = . Vj and Z denotes the "cross-reactive immune response directed against all 
different virus strains" (Nowak and May, 2000 p. 130). 

In the models they study, they set various parameters to values assumed to be constant. 
More precisely, they assume the parameters do not depend on N. Why is it impossible ? 
Assume there were only linear and non-linear terms like: 

^ = r NM Vj - p NM Tj Vj Vj. 

What is the integrated (in j) equation ? If we assume all the populations Tj and Vj are 
independent on j, then Tj = ^2 k Tk/N = T/N and Vj = V/N. Then the integrated (in j) 
equation is: 

dV T/ PNM rp t/ .nrp T T/ V 

= r NM V — T V lf T, = -,^ = -. 

Since it is the only equation biologically measurable, the measured parameter in an ex- 
periment will be Pnm/N and not Pnm- In other words Pnm = O(N). This should modify 
the mathematical study. 

The conclusion of their study is that in their models "the cross-reactive immune re- 
sponses provide a selection pressure against antigenic variation, while strain-specific re- 
sponses select for antigenic variation" (Nowak and May, 2000). 
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2.2 A reduced Snedecor's model 

In Snedecor (2003), the author gives three models of the multiplication of virus HIV-1. 
Her goal is to model both drug resistance, the different behaviors in the lymphatic tissue 
and the peripheral blood, but also the immune system. Distinguishing the blood and 
lymphatic tissues drives the author to have numerous constants that model the fluxes 
between these body parts and enriches her model. 

In order to have a common basis with other models, we need to reduce the above 
mentionned models to make their main features comparable with the features of other 
models. After non-dimensionnalizing the three fields with respect to the health value 
of lymphocytes (T* = 2.5 10 11 ), we are lead to only one model in which the uninfected 
lymphocytes are denoted by T, the infected ones by U and the virus by V. We take 
the values of the parameters in the same body part (lymphatic tissue). The reduced 
Snedecor's model is written: 

0(1 - T ) + ^+y( T -!)-(!- <*s)PsVT 
+(1 - a s )(3 s VT - aU (3) 
aU - a s V - (3 S VT. 

In this system, when the parameters used by S. Snedecor appear in terms that do not 
appear in other models, they are denoted with her notation with an index S. For instance, 
the division rate of T cells is r s = 0.004 day" 1 , the treatment efficacy is 0:5 G [0, 1], and 
the viral clearance is as = 2 day -1 . After non-dimensionalizing, the other parameters 
become 

Ps = 0.0125 day" 1 , 75 = 4 x 10" 5 ,/3 = 0.01 day -1 , a = 0.7 day" 1 , and a = 250 day" 1 . 




2.2.1 Fixed points 

In the present subsection, we prove the following proposition: 
Proposition 1. There exists a threshold 

a((3 s + G- s ) , 

as4 = 1 7r (4) 

aps 

such that if the efficacy parameter as is above 0:54 then health is the only fixed point. If 
as < CC54 there are two fixed points: health and a seropositivity. 

Proof. 

The search for fixed points gives two possibilities. 



The first one is health: T* = !,[/* = 0, V* = 0. 
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The existence of the second one depends on the therapy's efficacy parameter as- 
Different critical values of as will appear in the discussion. The solution for T is 

T * = <?s/!3s 

g(l-as) _ i ' 

a 

and is drawn in Figure [3] (left) with an infinite value for as = as3 = 1 — a/a = 0.9972. 
So, should a drug be very efficient (as > ct 5-3)7 then T* < and health would be the only 
solution. Then the solution V* is a non-negative solution of the second order equation: 

V 2 (l - a s )f3 s T* - V((r s - (3)(T* - 1) - (1 - a s )(3 sls T*) + f3 ls {T* - 1) = 0. (5) 

Numerically, it seems that the discriminant is non-negative (see Figure [3] right) but 
indeed, it is negative between 0:51 ~ 0.54920378 and as2 — 0.5492747378 and the min- 
imum is about —2 x 10~ 13 ! Notice that when perturbing the parameters this behavior 
remains. Except for as not too close to 1, the discriminant is small (see Figure [3] right) . 
We have drawn in Figure [3] (center) the only admissible solution V* as a function of as- 
It is then simple to have U* — (1 — a s )/3sV*T* /a. 




0.5 10 0.5 10 0.5 1 

Figure 3: The fixed points T* (lymphocytes), V* (virus) and the discriminant of ([5l) as a 
function of the treatment efficacy as- 

If we take care of retaining only admissible solutions ((T*, U* , V*) > 0), we must force 
T* > or as < «s 3 . Moreover, if T* crosses 1 (for as = «5 4 = 0.5492 exactly see ((!])), 
the sign of the constant term in ([5]) changes and so does one solution V of ([5]). One may 
then summarize the discussion for T*, and the two solutions V* and V£ of (j5J) in Table CD 
Indeed, there exists a non-health solution only for as < 054- Such a fixed point can be 
named seropositive solution. □ 

2.2.2 Stability 

Since the domain of admissible fixed points is not regular, the study of stability through 
the eigenvectors and eigenvalues at a corner (such as health) needs to be more precise. 
We define admissible directions to prohibit directions that do not enter the domain: 
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(T*,V*) 















Table 1: Second fixed point existence for reduced Snedecor's model. 



Definition 2. Let X i— > f(X) a smooth function and the associated dynamical system 
X'(t) = f(X). Let us denote B the biologically admissible domain (all biological fields non- 
negative). Let X* be a fixed point of a dynamical system (f(X*) = 0) at the boundary 
OB, and (A, u) one of its eigenvalue/ eigenvector. 

The eigenvector u is defined as admissible only if X > and either X* +eu or X* —eu 
for positive e enters the domain. 

The following proposition investigates the fixed points stability: 

Proposition 3. Let as4 as defined in (Jty. When health is the only fixed point (as > 
a si) it is stable. When there are two fixed points (as < asi), health is unstable in one 
admissible direction. 

Proof. The Jacobian matrix of the second member of ([3]) enables us to study the local 
behavior of the solutions. It is: 



J 



( r s V* 



f3-{l-as)(3 s V* -(1 - a s )(3 s T* + 



rslsiT* 



\ 



^ - OL S )PsV* 

-PsV* 



—a 
a 



(7s + V*) 
1 - a s )P s T* 
-as - (3 S T* 



(6) 



/ 



In the case of the first fixed point (health (T*,U*,V*) = (1,0,0)), the characteristic 
polynomial is —{(3 + A) (A 2 + \(a + crs + As) + a(o~s + Ps) — (1 ~ a s)o,l3s)- The discriminant 
of the second order polynomial is (a — cr s — Ps) 2 + 4(1 — as)aj3s > 0. So the roots are 
real and it is possible to discuss the sign of the roots. The already met threshold value 
«54 = 0.5492 (exactly see (j3J) is stil the key value for discussing on as- If as > 054, 
health is alone (cf. Proposition CD) and the roots are negative, so it is locally stable. If 
as < a s^ there are two stable and one unstable eigenvector. We easily check that the 
unstable eigenvector at this fixed point (health), located at a corner of the domain, enters 
the domain in the sense that one direction along this eigenvector (among the two) lets all 
the fields be non-negative. So this direction of instability is admissible. □ 



We checked numerically that the locally stable fixed point (as > 0154: health alone) 
remains stable even under non-small perturbations. In Figure HJ we draw the dynamics 
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of stable health. In this computation, we took the same parameters as Snedecor: 

r s = 4 x 10" 3 day -1 , a s = 2 day -1 ,/? s = 1.25 x 10" 2 day" 1 , 75 = 4 x 10" 5 , 
P = 0.01 day" 1 , a = 0.7 day" 1 , a = 250 day" 1 , 

and for the specific case of health stable, we took 0:5 = 0.6 > 0:54 for the treatment 
efficacy. The qualitative evolution of the pair lymphocytes- virus is comparable to the one 
of Nowak-May's model presented on Figure [H The number of virus loses a factor 5 in 3 
days typically. 




12 3 0.9992 0.9994 0.9996 0.9998 1 

time (days) uninfected lymphocytes 



Figure 4: Numerical results for Snedecor model with three equations. Case where the 
health is stable. Short time evolution on the left (3 days) and phase plane for lymphocytes 
and viruses on the right for a 600 days evolution. 

In the case where there exists also a second fixed point (seropositivity: as < 054), even 
odd initial conditions like (T , U , V ) = (1, 1, 1) lead to the second fixed point available. 
This second fixed point happens to be numerically locally stable as can be seen on Figure 
[5] where the parameters are the same as Snedecor's, as recalled above, except the therapy 
efficacy as = 0.3. In this figure, one may check that health is locally unstable as the 
solution evades health to get closer to a seropositivity. Notice that this second fixed point 
can be interpreted as seropositivity. But as it is stable, the model predicts no death ... 

2.2.3 Some comments 

The disappearance of the second fixed point (seropositivity) when as increases could be 
meaningful. Yet in real life, even for very efficient drugs, the virus kills. Despite highly 
active multitherapies, over ten percent of AIDS patients face therapeutic escape due to 
viruses which were indetectable for years while developing resistance against all approved 
drugs. So, to what extent can any stable fixed point be meaningful ? This criticism is 
valid for all the models analyzed in this article. 

The original model makes a distinction between the drug-resistant and drug-sensitive 
viruses which, in view of the goals of the author, is meaningful. But it makes no difference 
between infectious and non-infectious viruses. 
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50 100 150 200 250 . 4 . 5 0.6 . 7 0.8 0.9 1 

time (days) uninfected lymphocytes 



Figure 5: Numerical results for Snedecor model with three equations. Case where the 
seropositive state is stable. Relatively short time evolution on the left (260 days) and 
phase plane for lymphocytes and viruses on the right for a 900 days evolution. 

In a sense, the Snedecor's model takes into account the immune system. But it was 
proved that "the growth fraction of CD4 + (...) was correlated (...) with viral load" 
(Sachsenberg et ai, 1998). No such correlation appears in the model. Indeed, the immune 
system is modeled only through its exhaustion when V is too large. 

If as = (no treatment), the coefficient before the terms of disappearance of T, of 
V and appearance of U is the same. The identity of the V and T terms indicates that 
in this model, each time a virus infects a T, it disappears. So viruses are assumed to be 
free viruses. Indeed, the author defines her field V as "free virus". Similarly, the term 
—ftsVT in the evolution of V is by no means a model of the immune system response. 
Moreover, the author defines parameter (5s as "infection rate of ... T cells by ... virus". We 
emphasize here the coherence of Snedecor's model which takes into account "free virus": 
a free virus disappears each time it infects a lymphocyte. This property is not satisfied 
by most other models. 

In this model, viruses disappear only through natural death or infection of a T. So 
there is no immune system effect. Indeed, as is called "viral clearance (death) rate", but 
it does not depend on T. This is frequent in HIV modelling, but not realistic unless the 
immune system is neglected ! This is discussed in Perelson (2002) (pp. 31-32) where the 
author acknowledges his trouble: "The fact that models with constant [as] can account 
for the kinetics of acute HIV infection is surprising". 

2.3 Perelson's model 

In Perelson (2002) (and preceding papers with various coauthors), A. Perelson proposes 
a dynamical system to describe the interaction of HIV virus with T4 lymphocytes. The 
model uses four fields that we will denote with our notations to ease comparisons: un- 
infected lymphocytes (T), infected lymphocytes (U), infectious free viruses (V) and un- 
infectious free viruses (W). After non-dimensionnalizing with respect to the amount of 
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lymphocytes in a safe body (and a given volume), the system reads: 



XT = p(l-T)-5 P VT 

= SpVT-aU 
dj£ = aBU-apV 



(7) 



dW 

at 



a(l-6)U -a P W 



In this system, when the parameters used by Perelson appear in terms that do not 
appear in our model, they are denoted with his notation with an index P. 

2.3.1 Fixed points 

One may state the following theorem concerning the fixed points of . 

Theorem 4. There exists only two fixed points to system ffi). The first one is "health": 
(T*, U*, V*, W*) = (1, 0, 0, 0). The second one will be denoted seropositivity, where: 



T* = ^, V* = — (—-1 



W* = ^—^-V*, U* = ^-V*. (8) 
aO 



(9) 



The seropositivity is admissible (V > 0) under the condition that 

a5p9 — aop > 0. 
The proof of this theorem is easy and left to the reader. 

2.3.2 Stability of fixed points 

So as to evaluate the local stability of fixed points, one must compute the right hand 
side's Jacobian matrix J: 



(10) 



J(T*,U*,V*,W*) 
In the case of "health", (T*, U*, V*, W*) = (1, 0, 0, 0) it looks: 



/ 


-(3 - 5 P V* 





-5 P T* 


\ 




5 P V* 


—a 


5 P T* 










a9 


-Op 





\ 





a(l-0) 





-<T P J 



J(l, 0,0,0) 

It enables to state a theorem of conditional stability for "health". 



/ 


-13 





-5p 


o \ 







—a 


S P 










a9 


-a P 





V 





a(l-6) 





-Op J 



fill 
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Proposition 5. Health as a fixed point is stable if and only if 

a5 P 9-aap<0. (12) 

The proof is very simple and left to the reader. The eigenvalues are real and take the 
values —(3,—o-p, 1/2 (—a — op ± a/(« — o>) 2 + Aa5p9)). The admissibility of the unstable 
direction can also be checked. 

In the case seropositivity may occur, one may state the following theorem. 

Theorem 6. Under the admissibility assumption the seropositivity fixed point (M) is 
stable. 

Proof. Apart from A = — o>, the eigenvalues satisfy: 

A + + 5 P V* 5 P T* 
-5 P V* X + a -5 P T* 
-a9 A + o-p 



0. 



or A 3 + bi\ 2 + b 2 \ + h 
1 



with 



aop 



((3a5p6 + aap(a + o>)), 62 = — — —(a + ap), 63 = [3(a5p6 — aap). (13) 



OLOp 



We need to use the Routh-Hurwitz Criterion (Henrici, 1974 p. 490) that gives neces- 
sary and sufficient conditions ensuring that the roots of the cubic polynomial have positive 
real parts. In our case, this criterion reads: 



A x = bi > 0, A 2 



h 1 

h b 2 



>0, A 3 



61 1 
h b 2 b x 
63 



> 0. 



The first condition is obviously satisfied. Thanks to condition ([9l), the third condition 
is equivalent to the second one. It happens that A 2 can be computed: 



A, 



1 



a 2 ap 



(f3a5 P 8 + aa P (a + cr P ))(af35 P 6(a + a P ) + a ap) - f3(a5 P 8 - aa P ), (14) 



and this term is obviously positive because the only negative term is compensated by one 
of the expanded terms. □ 

2.3.3 Some numerical simulations 

So as to have simulations comparable with the other models studied in the present article, 
we use the same values of parameters as previously: 



P = 0.01 day" 1 , a = 0.7 day" 1 , a = 250 day -1 . 
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Moreover, some other parameters were the same as in the Snedecor's model and we used 
the same values: 

8 P = 0.0125 day" 1 , a P = 2 day" 1 . 

So as to be in the regime of health stable, we used 9 = 0.1. As can be checked on Figure El 
health is stable, although it needs numerous days to come back to health. The qualitative 
evolution of the lymphocytes and viruses is comparable to the one of Nowak-May's model 
presented in Figure Q] and to the one of Snedecor's Model in Figure 0J The uninfected 
lymphocytes decrease to their minimum value (0.999625) in about six days. Then some 
hundreds of days are needed to get sufficiently close to 1. 




Figure 6: Numerical results for Perelson's model with four equations. Case where the 
health is stable. Short time evolution on the left (3 days) and phase plane for lymphocytes 
and viruses on the right for a 600 days evolution. 

So as to have a stable seropositive state, we only change 6 for 9 = 0.6. Results can be seen 
on Figure [7j The global dynamics converges to the seropositivity with oscillations. The 
maximum value of the free virus (about 2.2) is reached at day 39. Then the uninfected 
lymphocytes decrease until 0.599 at day 58. The minimum value of virus (0.019525) 
happens on the 116th day and the cycle goes on as shown in Figure [71 



3 A new model with antigenic variable 

So as to build up our model, we will follow the biological description of the various 
phenomena concerning the various fields: the uninfected lymphocytes of antigenicity i G 
A (Ti(t)), the infected lymphocytes of antigenicity i (Ui(t)), the infectious viruses of 
antigenicity i (Vi(t)) and the non-infectious viruses of antigenicity i (Wi(t)). The space 
A is still undetermined. The best set is not known, but when trying to get a macroscopic 
model (by integration of antigenicity i), we will need to investigate the various possibilities 
to get a limiting operator for (t, i) G IR + x A. This is the issue of a forthcoming paper. 
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Figure 7: Numerical results for Perelson model with four equations. Case where the 
seropositive state is stable. Relatively short time evolution on the left (200 days) and 
phase plane for lymphocytes and viruses on the right for a 600 days evolution. 

We will also need the sum of each field: 

T = Y^T j ,U = Y i U j ,V = Y i V j and W = J^Wj. 

3 3 3 3 

3.1 Lymphocytes evolution 

The variation in time of Tj(t) (dTj/dt) must take into account various phenomena: 
• a natural death and generation modeled by a term like: 



dT) 
dt 



= -PjT s + 7, (15) 

natural 



• when a virus (Vj or Wj) of antigenicity j is detected, the immune system generates 
lymphocytes of the same antigenicity to fight them. This can be modeled by an expo- 
nential multiplication whose time constant is roughly proportional to the inverse of the 
number of viruses Vj + Wj. This is a "Lotka-Volterra" type term that we met only in 
Pastore and Zubelli (2008) for the cross-antigenic immune action: 

(^f) =+C 3 (V 3 + W 3 )T 3 (16) 

V Ut / growth 

• the virus attack lymphocytes independently of the antigenicity. More precisely, there 
are two regimes: 

+ If V/T < 1 as a virus may attack only one lymphocyte T k (Dern et al, 2001), the 
ratio of attacked lymphocytes is V/T. If we introduce a time constant, it gives a term 
like: 

J ' ^ T„ if - < 1. 



^ ) infection T 3 V T ) T 
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+ If V/T > 1, as a virus may infect only one lymphocyte (Dern et al, 2001) only part 

of the viruses may infect the lymphocytes. More precisely, no more than T viruses may 
infect. In other words, the ratio of attacked lymphocytes may not be superior to 1. There 
is a kind of saturation of the efficacy of predators Vj while the number of lymphocytes 
decreases with time. Up to a time constant, we saturate the infection ratio: 

V dt / infection T 3 * 

+ If V/T ~ 1, the two terms must match by continuity, and so Tj = fj. This complex 
behavior may be compiled due to a non-linear function J of the "min-mod" type: 

Mil- cm 

The effect discussed here can be modeled by 

V dt J infection T 3 \ T J 

In order to justify our essentially linear term in (TT8l) in an other way, we make hereafter 
various assumptions and wonder how our term should behave upon these assumptions, 
seen how reality behaves. 

Firstly we assume that we double Vj (and only it) without changing V or to a noticeable 
extent, while V/T remains small. Then the effect does not change and so the term must 
not change. This prohibits a simple term like —Vj/tj. 

Secondly we assume that all the VJ are doubled (including Vj) and so V is doubled 
too for a still small V/T. Obviously the effect is doubled. So the modeling term should 
depend on VJ (whatever i) only through V. 

Thirdly we assume that we double Tj, and only it (not the more general Tj), without 
changing T (or to a noticeable extent) nor V (and still V/T <C 1). Then the effect doubles. 
It proves that the term should depend linearly on Tj. 

At this stage, we have a linear dependance both in V and Tj. If we do not take care, 
we might deduce a —VTj/rj term that leads to —VT/t when integrated over j. We still 
need to explain why our term may not depend linearly on T although it depends on Tj. 

To that end, we make the assumption that we double all the T (including Tj). Since 
a virus may infect only one lymphocyte at a time, the effect should not be modified 
as the limiting parameter is not the amount of lymphocytes but the amount of viruses. 
Indeed, assume a patient has caught any disease that makes his immune system produce 
lymphocytes, it should not make the HIV infection more virulent, in the first stage where 
V/T <C 1. This is satisfied by our term and not by —VTj/rj. 

Notice that by summing over j and assuming Tj constant, we find dT/dt = — V/r in 
the first regime (V < T) and dT/dt = —T/r in the second regime (V >T). Moreover, 
if F > T the effect does not depend on V as the limiting factor is the presence of 
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lymphocytes. These are expected as they seem natural. Most authors have considered 
that the infection term should be quadratic like —VT/r. This conclusion may not be 
drawn from our considerations. It is even denied as the regimes exclude one another. To 
our knowledge, a term like ours has never been proposed. It must be recalled that all the 
models are designed to be tested with the only biologically meaningful fields T + U and 
V + W. Yet the phenomenon we model will simplify in the evolution equation of T + U 
and V + W. So it should never be measured. Yet, it translates into the model a crucial 
biological reality. So this term would deserve to be tested and such a test is postponed 
to a forthcoming research. 

The consolidated evolution equation for the lymphocytes is the sum 



dlj_ 
dt 

dT, 
dt 



dJ) 
dt 



+ 



natural 



d7} 
dt 



+ 



growth 



dJ) 
dt 



- i 3 jTj + lj + C j (V j + W j )T j --J 



infection 



(19) 



3.2 Infected lymphocytes evolution 

The evolution of infected lymphocytes depends on various effects: 

• the infection of a lymphocyte by a virus (same term as for Tj <HM with a plus sign) 
generates an infected lymphocyte Uf 



dU, 
dt 



generation 



1 ,(v 



and a natural death: 



To summarize, we have: 



dOj 
dt 



natural 



-ajUj. 



d^ 
dt 



dt 



+ 



generation 



d^ 
dt 



natural 



T 



(20) 



3.3 Infectious and non-infectious virus evolution 

The multiplication of virus is the main phenomenon and occurs in the infected lympho- 
cytes. So the sum in j growth term of V + W must be 

E p±M) (21) 

j V Ui / growth V ut / growth 

as taken into account by Nowak and May (2000), Snedecor (2003) and Perelson (2002) 
(whether they distinguish infectious and non-infectious virus or not). 
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The evolution of virus depends also on mutations. Some of them modify the abil- 
ity to infect and others the antigenicity. Here, those two properties are considered as 
independent, as mentionned in the introduction where references are given. 

Only one parameter 9 measures the probability to mutate to an infectious offspring 
(necessarily from an infectious virus). So 1 — 9 is the probability to mutate to a non- 
infectious offspring. Since we assume antigenicity and ability to infect are independent, 
we may assume 9 does not depend on j (nor k). 

Let us denote Skj the probability to mutate from an antigenicity k (only Vk as Wk 
does not even infect and so does not mutate) to an antigenicity j (either Vj or Wj) per 
unit of time. So we have: 



and Skj depends a priori on the number of antigenicities. The offspring of a virus Vk will 
mutate to Vj with the probability OSkj (the case k = j where there is no mutation is 
included). So it will mutate to Wj with the probability (1 — 9) Skj- 

As a consequence of the assumption that antigenicity and ability to infect are inde- 
pendent, equation (T2TT) can be split into two parts : 



Notice that the value of 9 suggested by Thomas et al. (2007) ranges from 1 to 1/60,000. 
But Kothe et al. (2006) found in one experiment a ratio 9 close to 1/8. To fix the ideas, 
we suggest in the following to set 9 = 1/10. 

The value of Skk is considered as relatively low by the biological community but we have 
not found any precise and trust-worthy proposition in the literature. 
The distribution k \— > Skj for j ^ k needs to be non-uniform for a biological reason. As is 
well known (Bebenek et al, 1993; Ji and Loeb, 1994), there exist mutationnal hot spots 
in the HIV virus. We may guess these hot spots will elicit non-uniform distribution of the 
antigenicity during multiplication/mutation. The precise value of the probability k \— > Skj 
is of course an open problem. 

3.3.1 Infectious virus 

The infectious virus variation (dVj/dt) depends on various effects. 

• The first and most complex is multiplication and mutation. Through mental experi- 
ments, we are to determine the fields present in the modeling term. 

As multiplication takes place in the infected lymphocytes, if there were no such lym- 
phocytes (U = 0), whatever might be the number of (free) viruses, there would be no 
multiplication and the term would be zero. This would almost occur in the end of the 
disease. More precisely, as once a virus has infected a lymphocyte, no more viruses may 





(22) 




(23) 
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infect it (Dern et al, 2001), U is a good measure of the total number of infectious viruses. 
Note that in our model V is the number of free and infecting viruses. So there must be a 
dependance on U (or Uj at this stage). 

As the antigenicity of an infected lymphocyte has no link with the antigenicity of the 
multiplicating virus, the involved field must be U (and not Uj). 

As there are mutations from any antigenicity k, the Vk mutate and their offspring is 
made of either Vj or Wj with probability OSkj or (1 — 6)Skj respectively. So there must 
be J2 k SkjOVk in the modeling term. 

What can be the modeling term ? Up to now, we have U(J2k Skj@Vk) times an 
unknown term A(t). Once summing over the antigenicity, we must find aOU . So because 
of J23J), A(t) is such that: 



E (ir) = E A ® u ( E s * ev >> ) = a6U =► A & 

o V / growth j \ k I 



a 



because ^ ■ Skj = 1. So our modeling term for multiplication with mutation to Vj (V* 

dVj\ u 



Vj) is 



df , a-J2S kj 0V k . (24) 

/ growth k 

attack of the viruses by the lymphocytes of the same antigenicity produced by the 
immune system: 

/<U '^ -SjVft. (25) 



/ lymphocyte 

Once compiled, the evolution equation is: 

V / growth V / lymphocyte k 

3.3.2 Non-infectious viruses 

The variation of non-infectious viruses (dWj/dt) may be modeled by various effects and 
translated into various terms similar to infectious ones: 

• mutation from an antigeniticy k (Vk with Vj included) to Wj-. 



growth k 

• attack of the virus by the lymphocytes of the very same antigenicity. Notice that as 

the immune system may not detect whether a virus is infectious or not the £ 3 - must be 
the same as for Vj\ 



dt 



= -t j W j T j . 

lymphocyte 
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Once compiled the evolution equation reads: 



\ / growth V / lymphocyte y. 

As for Vj, the lymphocytes do attack only the virus of the same antigenicity. Moreover, 
they cannot make a distinction between infectious or non-infectious viruses. 



3.4 The dynamical system 

When we collect equations (TT9l) . (1201) . (1201) and (12T1) . our model reads finally as said in 
Dubois et d. (2007): 



<l/j -/3 J T,+ 7j + Q(^ + ^)T,-ij^)T„ (28) 



dt 7 ' J " J J/ J r/ V'/' 

^ - ^(D^-aM, (29) 



dt r, V T 

dVj aOU 



dt V 

k 



^.StiVt-dVjTj, (30) 



fc 

Very simple manipulation enable macroscopic laws to be proven: 



(31) 



d{T + u + Y J A^-iy 3 + w J ))) _ _ 



E + E Ti - E + af/ E ir- ( 32 ) 



dt ^- ; ' J ' ^ ,J ^- ; ' ' ~" ^ C, 

3 3 3 3 

Since we see no reason why ^ or Cj should depend on j, the equation (1321 ) could be 
considered as a simple linear combination of the integrated versions of (128H31I) . Such a 
law could be experimentally checked. 

There remains the problem of initial conditions and the way they enter into the evo- 
lution model with mutations. Indeed, in most ordinary differential systems, as soon as 
a function is identically zero in a subdomain, it remains so (Cauchy-Picard Theorem). 
Moreover, 10~ 20 or 10~ 40 are numerically very different while they both mean "zero". In 
simulations of mutation, we will need a quantic jump. All this is postponed to a forth- 
coming article. 
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4 Mathematical properties 



Starting from ( I281I3T1) . we make the various fields dimensionless by using the value of 
^equii = a t equilibrium as a characteristic value. Then we let N — 1: 

Ta Uj Vj Wj 



rji ' rji ' rji 7 rji 

equil equil equil equil 

We also define dimensionless parameters: 

^ Cj Tequih C -^equil- 

With these notations, the system reads: 



^ = P(l-T(t))-^J (J^j +u;(V(t) + W(t))T(t), (33) 
dU T(t) fV(t)\ 

^ = a^(t)-C^)^), (35) 
dW 



dt 



a{l-6)U{t)-(W{t)T(t). (36) 



The unknown fields for ( 13311361 1 are dimensionless whereas the time variable remains di- 
mensioned (by days). Since we assume there is only one antigenicity (N = 1), mutation 
disappears. Hereafter, we prove the solution remains admissible and exists globally. Then 
we look for the fixed points and study their stability. 

4.1 The solution remains admissible 



We intend to prove that the system (I33H36I) is mathematically well posed and biologi- 
cally meaningful: it does not exhibit negative values of the fields for admissible initial 
conditions. To that purpose, we make various assumptions on the parameters: 

P, r, u, a, ( are real positive and < 9 < 1, 

J(») is a real function concave over [0, oof, J(0) = 0, J'(0) = 1, (37) 
J(x) — > 1 when x — > +oo and J is bounded on Mr. 

>From the biological meaning of the system, we define the set of admissible fields: 

B = {(T, U, V, W), T > 0, U > 0, V > 0, W > 0}, (38) 

and its interior: 

B = {(T, U, V, W), T > 0, U > 0, V > 0, W > 0}. (39) 
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Due to the Cauchy-Picard theorem, we know that there exists a local in time solution. 
The question is then whether this local solution is admissible. Our main result is the 
following Theorem. 

Theorem 7 (Biological consistance). The solution (T(t), U(t), V(t), W{t)) of system (E 
E2IJ with an initial condition in B remains in B. 



To prove Theorem 0, we will need various lemmas. The first one states that the number 
of lymphocytes does not vanish in finite time. 

Lemma 8. If the initial condition (T , Uq, Vq, Wq) is in B, then T(t) > for all time for 
which the fields T, U, V, W are defined. 

Proof. Thanks to the Cauchy-Picard theorem (Arnold, 1978), we know that for t suffi- 
ciently small, T(t) is positive. So if there exists at least one time t such that T(t) = 0, 
then we define t* to be the smallest and we have t* > 0. Since J(«) is bounded (whatever 
V and T), (1331 ) writes at that time 

f(O=/»0. 

As a consequence, for t below and sufficiently close to t*, T(t) < 0. But as T > 0, from 
the intermediate value theorem, there exists a t' smaller than t* for which T(t') = 0. This 
contradicts the assumption that t* is the smallest and completes the proof. □ 

We need now to study the various cases where the initial conditions are either in B or 
on the boundary of B. This will be discussed through some lemmas where the initial 
condition has either zero, three, two or one initial vanishing fields. In the case of initial 
condition in the interior of B, one may state the following Lemma. 

Lemma 9 (Zero vanishing initial condition) . If the initial condition of system [SWSh}) is 
in B, then the solution remains in B for any t > provided it exists. 

Proof. We will discuss the cases where one, two or three fields vanish simultaneously. 

• Let us assume W vanishes first and alone (before the other fields) and then let us 
denote t* the smallest time for which W vanishes. On [0, t*], one has: 

V(t) > 0, U(t) > 0, W{t) > 0, W[t*) = 0, (40) 



AW 
At 



in addition to the fact that T(t) > (see Lemma ED- The equation ( 1361 ) writes 
a (1 — 9) U(t*) > because < 9 < 1 (see the assumption (1371)). So, there exists a time 
t smaller than t* at which W(t) < 0. From the intermediate value theorem, one may 
conclude that there exists a time smaller than t* (and than t) at which W vanishes. This 
contradicts the definition of t* and so W may not vanish first. 

Identical arguments enable to prove that V may not vanish first. To prove that U may 
not vanish first neither, we assume that t* is the first vanishing time of U. So, on [0, t*]: 
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U{t) > 0, U{t*) = 0, V(t) > 0, W{t) > 0. (41) 

Thanks to LemmaEJ one has T(t*) > and so ^(t*) = J(^) ^ > 0. In a similar way 
to the two previous cases, one gets a time t smaller than t* for which U(t) < 0. Thanks to 
the intermediate value Theorem, one gets also a time smaller than t* where U vanishes. 
This contradicts the definition of t*. So U may not vanish first neither. 

• Let us discuss now the case where two fields vanish simultaneously. The case where 
V and W vanish simultaneously and alone is impossible for the same arguments as the 
case where W vanishes first. The case where U and W vanish simultaneously and alone 
may be treated in the same way as the case of U vanishing first. The only remaining 
case is if U and V vanish simultaneously (and not W). Let us then define t* the smallest 
such vanishing time. The Cauchy-Picard Theorem, applied to the system f!33H36f) with 
reversed time, enables us to claim that the (unique) solution is also such that 

dT dW 

— =P(l-T(t))+uW(t)T(t), JL = -{ W (t)T{t), U(t) = V(t) = 0, 

for all time t in [t* — e, t*] for some e > 0. This contradicts the definition of t* as the 
smallest vanishing time. 

• In the case where U, V and W vanish at the same time, the proof is the same as for 
Lemma [T0l □ 

The following Lemma solves the case of three vanishing initial fields. 

Lemma 10 (Three vanishing initial conditions). If the initial condition is T > 0, Uq — 
Vq = Wo = 0, then the solution is unique and is health. 

Proof. The proof relies only on the Cauchy-Picard Theorem which states that 

T (t) = 1-(1- T ) exp (-/%), 

and the other fields identically zero. □ 

The next Lemma deals with the case where two initial fields vanish. 

Lemma 11 (Two vanishing initial conditions). If the initial condition is among 

T > 0, U = V = 0, W > 0, (42) 
T > 0, U > 0, V = Wo = 0, (43) 
T > 0, U Q = 0, V > 0, Wo = 0, (44) 



the solution of system ftEMEE) remains in B. 
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Proof. 

• In the case (l4*2l) . the Cauchy-Picard Theorem enables us to claim that there exists a 
solution of (EMSD for which U(t) = V(t) = 0, and the fields T(t), W(t) satisfy 

^ = f3(l-T(t))+uW(t)T(t) (45) 
dW 

— = -(W(t)T(t). (46) 

Let us assume W vanishes at some times. Among these times, we chose t* to be the 
smallest. On the compact set [0, t*], T is continue and so bounded by 7 > 0. So f j46"l ) 
enables to claim 

dW 

> then W(t) > W exp (-(yt) for t > 

(J.6 

which contradicts the assumption that W{t*) = 0. So W may not vanish in finite time 
and so in the case ( 1421) . PU(t) > for any positive time and U — V = 0. 

• In the case (iHD, % > and so V(t) > for t small enough. Similarly, for t small 
enough, W^(t) > 0. So for t small enough, the solution enters the interior of B. Such a 
new initial condition has been adressed in Lemma EE So the solution remains in B. 

• In the case (HH), > thanks to (EH), and so U(t) > for t sufficiently small. 
Indeed, ^ = but ^ = a (1 - 9) ^ > 0. This is enough to assess that W{t) > for 
t sufficiently small. Using the new "initial" condition at this time, we are driven back to 
the case treated by Lemma [9l This completes the proof. □ 

Lemma 12 (One vanishing initial condition). // one and only one initial field vanishes, 
the solution remains in B. 

Proof. If U = (and V W > 0), then ^(0) = J(^) > 0. One may conclude in 
a similar way to the case (1441 ) treated in Lemma [TU If Vq — (and Uq Wq > 0), then 
tt-(O) = a9Uo > 0. So in finite time, one is driven back to the case treated in Lemma 
M (zero vanishing intitial condition). If W = (and U V > 0), the argument is very 
similar. □ 

Proof of Theorem [3 

Up to now, we have proved that if the initial condition does not vanish (Lemma I), 
vanishes three times (Lemma QUI), two times (Lemma [□]), or once (Lemma d2|), the 
solution remains in B. This completes the proof of Theorem [71 □ 



4.2 Global existence 



The following Theorem states that the solution remains finite and so is global in t. 
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Theorem 13. Let 

rj = CLLO-aC, (47) 

and 



7 



(48) 



// the initial condition (T , Uq, Vq, Wo) is in B, then the solution of f3ffl3b)) satisfies: 

T(t) + U(t) + ~(V(t) + W(t)) <T + U + ^(V + W ) + pt if V <0, (49) 
and 

T(t) + U(t) + ^(V(t) + W(t)) < (T + U + ^(Vo + W ))e^ + ^(e^-l) if V > 0. (50) 

As a consequence, the solution is finite for all time t G [0, +oo[ and so is global. 

Proof. By adding (13311341) and j times the sum of (l35l) and (l36l) . all the non-linear terms 
disappear and one has: 

Thanks to Theorem [2, the solution remains in B and so U > 0. Moreover if rj < 0, we 
can minorate the fourth term of (ED) by because ( > 0. As a consequence, 



A(V + [/ + ^(V + M0) < /3, 

and the relation (T49l is a simple consequence of the integration in time of the previous 
inequality. 

If rj > 0, we integrate between and t the equation (I511 . and owing to the positivity 
of T, one has 

„t 

T + U+^(V + W)<T + U + ^(V + W )+pt + 1 j U(t')dt'. (52) 



o 



Let us denote 



5 = T + U + j(V + W ) and <j>(t) = / C/(t') dt'. 

Jo 

Due to Theorem [7] we know that the solution remains in B (T > 0, V > 0, W > 0), the 
equation ( 1521 enables to write: 

^ < 5 o + ^ + 70(*), 
A(exp(- 7 t)0) < exp(- 7 t)(5 + /3t). (53) 
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The inequality (1531) can be integrated (0(0) = 0): 
exp (--yt)<j>(t) < \--t - - ( 5o + - 



7 7V 7 



x 1 ft P 

exp (-7t) + - [5 + - 
7 V 7 



or 

7 0(t) < - (pt + 6 + + \ S + exp (+ 7 t). 
With such a bound for 0, one may take back the right hand side of (io2l) : 

T + [/ + %{V + W) < 5 + pt + 70 < S exp ( 7 t) + 2 (exp ( 7 t) - 1), 
C 7 

which is the inequality (l50l . This completes the proof. □ 

4.3 Fixed points 

Looking for fixed points of ( |33ll36l) . we must find solutions of the associated stationnary 
system: 

= p(l-T)-£jn£\+u(V + W)T, 

T J (t) = aU > (54) 
adU = (VT, 
a(l-9)U = (WT. 

We will prove the following proposition 
Proposition 14. Let rj defined in IJ1\) and 

auj — a^' 

The fixed points of {Sffliffl) depend on the sign of n. Three cases must be distinguished. 

• lfr)>0, the fixed points are either health and seropositivity (if p < 1), or only health 
(ifp>l). 

• Ifrj — 0, there is no fixed point else than health. 

• If r] < 0, health is always a solution. Moreover, there appears a non-explicit threshold 
value L and three sub-cases depending on L: 

+ When p < 1 there is one seropositivity fixed point. 

+ When 1 < p < L, there are two seropositivity fixed points. 

+ When p > L there is no seropositivity solution. 
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Proof. One finds easily that 



(1 - Q)v = ew, 
U=iiV + W)T = -^VT. 



(56) 



Then, the system (IBTl) reduces to 





(3(1 - T) - ^VT + ^VT, 



a9 



VT. 



(57) 
(58) 



Since T = is not a solution, one may simplify T in (15811 . Equation (1571) gives T: 

(3 1 



T 



Thanks to the value of T given by ( l59l) . we are driven to solve (1551 ) in the form: 



(59) 



7 



a0 



■V = pV. 



(60) 



Three cases appear to solve this equation. 

• If n > 0, we must find the solution V of (l60l) where V > 0. To guarantee T > 0, we 
must have V < V. Since the function V i— > J(V^(V — V)/V) is symmetric with respect 
to V/2 and the right hand side pV is a linear function of V, only two subcases must be 
distinguished according to the Figure El 




v v 



Figure 8: Two characteristic shapes of the curves V 
(dotted) for rj > 0. 

This discussion is summarized in the following: 
1. If p < 1 there are two solutions: 

(i) Vl = 0^T 1 = l,U 1 = = W 1 ; 

(it) V 2 >0^T 2 >1. 



pV and V ^ J(V(V - V)/V) 
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The first one will be denoted health and the second one seropositivity. We notice that 
V2 < V and in the end, the condition on V to ensure T > is satisfied. 
2. If p > 1, the only solution is health. 



• If n = 0, then V = 00. There is no solution. 

• If rj < 0, then V < and there is no constraint on V to ensure that T > 0. Unlike the 
case r] > 0, the parabola inside the function J is of the type y = +x 2 instead of y = —x 2 . 
So as to circumvent this, we will invert the function V 1— > V(l — V/V) on M + where it 
is one-to-one. For any V G M + (we look for non-negative V), there is a unique X G M + 
such that 

- V '+ \ V 2 -AVX 
X = V(l - V/V) V = ip(X) = L . (6i) 

Then equation (l60l) amounts to 



J(X) = (4{X). (62) 
In other words we need to intersect J(X) and the parabola pip(X). Only three cases 




Figure 9: Three characteristic shapes of the curves X \— > pij)(X) and J (dotted) with J = 
tanh and n < 0. 

must be distinguished, according to Figure El They depend entirely on the slope of fnjj{X) 
at X = whose value is p. Hereafter, we denote with a superscript the roots of ( |62l) . 

1. If p < 1, there are two solutions denoted X}: 

(i) X\ = corresponding to health ; 

(ii) X\ > corresponding to seropositivity. 

2. If 1 < p < L there are three solutions denoted by Xf: 

(i) X\ = corresponding to health ; 

(ii) Xf > corresponding to seropositivity ; 

(Hi) < Xf < Xf corresponding to seropositivity. 
Notice that there is no explicit value of L. Yet for J = tanh and V ~ —0.054, we could 
numerically find L ~ 3.7. 

3. If p > L, there is only one solution Xf = (health). □ 
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4.4 Stability 

The Jacobian matrix is, at every fixed point: 
/ 



Tt It 



T 



\ 



J (t)/ t Tt J ' (t 
-CV 

-<;w 





—a 
aO 



a(l 



(jT-J'Qpj/r uoT ^ 

J' (£) It 

-CT 

-CT J 



If we compute the characteristic polynomial through the last right column, we find: 



P(A) = -(CT + A) 



-uT 



J 



J 



V_j> (V 

Tt \T 

V s 



/) ' Tt 

-CV 
-CW 



T 

VjffV 







ujT-J 1 



J' 



V 



T 



—a — A 

a6 -CT - A 

a - A J'(V)/t 



aO 
a{l-6) 



-CT-A 




(63) 



Thanks to (l56l . we can easily simplify the second line of the second determinant of (l63l) . 
As a consequence one may write: 



P(A) = -(CT+A) 



Y_ v (V 

Tt 
V 
T 



l_ x 

T J T 



Jil)/r 



V j,(V 



-ljT((T + \) 



J 



Tt 
(V 



T 



) I* - (*) 



o lut-j' np\ i 

a - A J 1 (f\ J j 
a8 -CT - A 

—a — A 
a(l-9) 



T 



In the general case, no more factorization could be found and the polynomial is 



P(A) = -(CT + A) 



\Tt 



J' 



I 
T 



-A 



+ 



f^-J't^/r) (a6(j(V)/T 



V 9 



T 



a6J' (}p\ It + (a + A)((T + A)) + 
-(a + \)(v]~ 



V-ji(V 
Tt It 



(64) 



4.4.1 Health's stability 

In the case of health (T = 1,U = = V = W), the eigenvalues are zeros of (1641) which 
can be rewritten thanks to (fool) : 



P(A) = (CT + A) 09 + A) ( ^(-1 + p) + A(a + C) + A ; 



0. 



(65) 



We will prove the following theorem. 
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Theorem 15. If p > 1, health is (locally) stable. 

If p < I, there exists one positive eigenvalue associated to an admissible eigenvector in 
the sense of Definition® 

Proof. The discriminant of the non-reduced second order polynomial in ([65]) is (a — () 2 + 
Aa6/r > 0. This enables us to claim that: 

• if p > 1, the four roots are negative and so the evolution is (locally) stable; 

• if p < 1, there exists one (and only one) positive root and so the evolution is locally 
unstable in one direction. 

We still have to prove that the eigenvector v associated to a positive eigenvalue of ( 1651 ) 
is admissible. In other words, we need to prove that v is such that for e > or e < 
small enough, (1, 0, 0, 0) + ev has its four components non-negative. To prove this, we 
write the system satisfied by the eigenvector: 



/ 


-(/3 + A) 





uj - 1/r 


u \ 


/ Xl \ 











-(a + A) 


1/t 

















a6 


-(C + A) 












\ 





o(l -9) 





-(C + A) / 


\x 4 J 




V / 



(66) 



where p < 1 and A is the (unique) positive root of ( |65i l. More precisely, A is the unique 
root of the third term: ad{— 1 + p)/r + \(a + () + A 2 . As this second order polynomial is 
the determinant of the 2x2 submatrix in the center of the matrix in (l66l) . and because 
of the very particular shape of the lines 2 and 3, we can claim these lines are bound. It 
suffices then to take off the third line to be driven to the system equivalent to (l66l) : 

(P + X)xi = (uj - 1/t)x 3 + LUX 4 
(a + X)x 2 = x 3 /t 
a(l-6)x 2 = (C + A)x 4 . 

As A > 0, we can see that there exists solutions such that x 2 ,x 3 ,x 4 be non-negative. So, 
at least locally the solution in the direction of this eigenvector is admissible and the proof 
is complete. □ 

4.4.2 Seropositivity's stability 

We have no rigorous study of the stability/unstability of seropositivity. So we use numer- 
ical simulations. In this subsection, we take 

P = 0.01 day" 1 , a = 0.7 day" 1 , uj = 0.01 day" 1 , a = 250 day" 1 , 9 = 0.1, 

and initial values 

T(0) = 1, U(0) = 0, V(0) = W(0) = 0.05. 

The other parameters (r, () are taken so as to illustrate the fixed points depicted in 
Proposition [HI We provide the evolution on a short time and a phase portrait for each 
case. Let us remind the fixed points: 
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1. If ?7 > 

(i) If p < 1: health (partially unstable) and one seropositivity numerically stable 
as can be seen on Figure [TOl In this simulation, n = 1.8, p = 0.28. The effect of the 
infection is to emphasize the activity of the immune system ; 

(ii) If p > 1: only health (stable), as in Figure [TT] (n = 0.4, p = 1.68). Minimum 
value for lymphocytes is obtained for t ~ 0.8 day and maximal one for t ~ 16.4 days. 

2. If ?7 < 0: 

(i) If p < 1: health (partially unstable) and one seropositivity numerically stable 
as can be seen on Figure IT2l (77 = — 4.5,p = 0.28). The effect of the infection is to reduce 
the activity of the immune system at a very low level (~ 0.015 times the level of health); 

(ii) If 1 < p < L: health (stable) and two seropositivities. We choosed the 
parameters ( = 6,r = 6 (77 = — 4.5, p = 1.008) for simulation. A first seropositive state 
(T* ~ 0.129, U* ~ 0.03073, V* ~ 0.992, W* ~ 8.9) was numerically found to be locally 
stable (not shown). A second seropositive state (T* ~ 0.992, f7* ~ 0.00028, V* ~ 
0.00117, W* ~ 0.01058) is very close to health and is locally unstable since it has three 
negative and one positive eigenvalues. So as to illustrate this, we have taken initial data 

X j = X* + eV j , 1 < i < 4, 

where X* is the fixed point, Vj is one of its eigenvector and e is sufficiently small. This 
experiment is depicted in Figure |T3l Note that although health is locally stable, initial 
conditions with a viral load of 5% drove the state to a seropositivity state. This proves 
that the basin of attraction is small. So we provide a simulation with an initial viral load 
of only 1% in Figure [H] ; 

(in) If L < p: health is stable as one may see on Figure [15] (77 = — 4.5, p = 2.8). 
Minimum value for lymphocytes is obtained for t ~ 2.6 days. 




Figure 10: Numerical results for our model. Case where a seropositive state is stable 
(r = 10, ( = 1). Relatively short time evolution on the left (25 days) and phase plane for 
lymphocytes and viruses on the right for a 600 days evolution. 



A Mathematical Model for HIV Dynamics 



33 




Figure 11: Numerical results for our model. Case where the health is stable (r = 20, ( = 
3) . Short time evolution on the left (3 days) and phase plane for lymphocytes and viruses 
on the right for a 600 days evolution. 




Figure 12: Numerical results for our model. Case where a seropositive state is stable 
(r = 1, ( = 10). Relatively short time evolution on the left (25 days) and phase plane for 
lymphocytes and viruses on the right for a 600 days evolution. The effect of the infection 
is to reduce drastically the efficacy of the immune system. 
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Figure 13: Numerical results for our model with 4 equations (r = 6, ( = 6). Simulations 
are initialized with states very close to the fixed point in directions that correspond to 
eigenvectors. The seropositive state is unstable in one direction among four. 
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Figure 14: Numerical results for our model with 4 equations. Case where health is stable 
(r = 6, ( = 6) only for small perturbations. Relatively short time evolution on the left (3 
days) and phase plane for lymphocytes and viruses on the right for a 1500 days evolution. 
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Figure 15: Numerical results for our model. Case where the health is stable (r = 10, ( = 
10). Short time evolution on the left (3 days) and phase plane for lymphocytes and viruses 
on the right for a 600 days evolution. 

4.5 Some comments 

Various effects are supposed to be more or less incorporated in any model and specifically 
ours: 

• If a model considers the field of free viruses, then infection makes an uninfected 
lymphocyte and a free virus disappear and an infected lymphocyte appear at the same 
time. Three identical terms (up to a ± sign) should be present in such a model. It is not 
often the case. Since we consider not only free viruses but also infecting viruses, we are 
not forced to have the same three terms. 

By explicitly deriving the model, we justifiy our "piecewise linear" term, although all 
the authors use a quadratic term. Very likely, the reason why this phenomenon's modeling 
has not been much studied is that the only biologically measured field is T + U. So the 
term modeling this phenomenon disappears in any evolution equation on T + U. Yet it 
models a crucial reality and deserves more attention. 

• We characterize each virus by its antigenicity and by the information that it is in- 
fectious or not. We characterize the lymphocytes by the virus antigenicity against which 
they have been designed. Mutation is then only a probabilistic phenomenon and the main 
modeling question is the space in which it takes place and its probabilistic law. Such a 
study is postponed to a forthcoming article. 

• We model the production of Tj by the immune system once it has detected the Vj. 
The production depends on the Vj population, and so our term is quadratic. It has no 
counterpart in any model we read, but it could be balanced by the quadratic term modeling 
the immune effect against each strain of virus present in numerous models (Pastore and 
Zubelli, 2008; ...) including ours. 

• The multiplication/mutation is the most challenging phenomenon. The death rate 
of infected lymphocytes is about 70 times greater than that of uninfected lymphocytes. 
Indeed, in Snedecor (2003), the author proposes various articles among which Finzi and 
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Siliciano (1998) for the death rate of infected lymphocytes and Sachsenberg et al. (1998) 
for the death rate of uninfected lymphocytes. The ratio of these rates is about 70. More- 
over, since only one virus may infect a lymphocyte (Dern et al., 2001), we may assume 
U is a good measure of the total number of infecting viruses. This assumption is very 
likely but would deserve to be further tested. A simple modeling of mutation generates 
our term. 

• We model the effect of the immune system against the viruses by a term depending 
quadratically on Vj and 7} since these terms are effective in the same regime. Such a term 
can be found in the models depending on antigenic variation in Nowak and May (2000) 
(chapter 12 and 13) and Pastore and Zubelli (2008) but in these models the lymphocyte's 
generation is modeled only through a linear term in Vj. Notice that with further assump- 
tions, one may find a linear combination of the T + U and V + W evolution such that 
this new combination is simply linear. This could be experimentally tested. 

The overall behavior of all the systems studied above (including ours if N — 1) al- 
lows the fields to remain non-negative and be attracted by some fixed points. So all 
these models predict convergence to some fixed point which is never immuno-suppresed 
(T = 0). We consider this to be a major drawback for the long term modeling of HIV 
infection. This opinion is shared by the authors of Pastore and Zubelli (2008) and they 
propose modifications to the Nowak-Baugham models enabling the longer term evolution 
modeling. 

What is the advantage of our model ? 

Notice that the only physical field is T = J2jTj and not 7\ (when N = 1). As a 
consequence, the widening of antigenicity support is a phenomenon not included in the 
case N = 1 nor in any other well-known "macroscopic" model reviewed above. 

Since our model takes the antigenicity into account, its macroscopic version should be 
more realistic. 

5 Conclusion 

We have thoroughly studied different models of HIV multiplication by systems of differen- 
tial equations. Some of them were reduced to be single-antigenic. With such a reduction, 
all of these models have fixed points that prohibit modeling of the last phase of the disease 
where the T count vanishes. This is also criticized in recent research (Pastore and Zubelli, 
2008). 

Moreover, we propose in this contribution a new model taking into account new phe- 
nomena among which lymphocytes generation by the immune system according to the 
presence of specific viruses, and immune effect against each virus strain. We model infec- 
tion and mutation/generation through new algebraic terms. This model is derived due to 
explicit arguments. It will be tested further in forthcoming research. 

Although the reduced version of our model has the same drawback of not enabling the 
immunity exhaustion, its general version takes into account the strains' diversity (here 
denoted as antigenicity) and the specificity of the immune response. So our full model can 



A Mathematical Model for HIV Dynamics 



37 



account for the last phase of the HIV infection where the lymphocytes' count vanishes. 
This will be studied in a forthcoming article. 
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